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INTRODUCTION 

This  project  addressed  the  FY10  PRMRP 
subject  of  epilepsy.  Increasing  incidence  of 
traumatic  brain  injury  (TBI)  among  soldiers 
will  likely  lead  to  elevated  levels  of  disability 
due  to  TBI-related  seizures  and  epilepsy.  The 
lack  of  a  reliable  biomarker  hinders  efforts  to 
interrupt  the  evolution  of  epilepsy  from  TBI. 

A  new,  network  paradigm  for  analysis  of 
brain  imaging  data  suggests  a  new  direction 
for  diagnosing  brain  injury.  Existing  analyses 
have  neither  been  applied  to  epilepsy  nor 
have  been  validated  by  gold  standard  data. 

The  purpose  of  this  research  was  to  initiate 
exploration  of  the  concept  that  network 
properties  of  imaging  data  within  the  scope  of 
predicting  the  transition  of  TBI  to  epilepsy. 

The  specific  objectives  of  this  project  were  to 
develop  a  fast  analysis  protocol  and  validate 
the  analysis  with  gold-standard  invasive 
electrophysiology  measurements  from 
epilepsy  patients.  The  innovative  aspects  of 
the  research  are  application  of  a  partial 
differential  equation  framework  for  fast 

analysis,  the  application  of  the  network  paradigm  to  epilepsy  patients  and  validation  with  gold  standard  invasive  measurements.  The 
relevance  of  the  project  to  the  FY10  PRMRP  topic  of  epilepsy  stems  from  the  potential  use  of  a  validated,  network  paradigm  as  a 
biomarker  of  risk  for  the  transition  from  TBI  to  epilepsy  with  the  intent  of  interrupting  that  transition. 
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Figure  1.  Failure 
of  deterministic 
tracking  for 
application  to 
intracranial 
electrodes. 
Deterministic 
tracks  (pale 
yellow  lines, 
indicated  by 
arrows) 
originating  in 
cortical  regions 
(A,  B)  fail  to 
connect  the 
cortical  regions. 
PDE-generated 
track  (red  line) 
readily  delineates 
connection. 
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Research  accomplishments  are  summarized  below  along  the  lines  of  The  Statement  of  Work,  which  outlined  the  following  main  tasks: 

•  Develop  a  partial  differential  equation  (PDE)-based  tractography  methodology  to  enable  fast,  whole-brain  measurements  of 
connectivity. 

•  Validate  noninvasive  measurements  of  connectivity  by  comparison  to  gold  standard,  invasive  electrophysiology 
measurements. 

•  Summarize  and  publish  results. 

Develop  a  partial  differential  equation  (PDE)-based  tractography  methodology  to  enable  fast,  whole-brain  measurements  of 
connectivity. 

A  principled  theoretical  framework  for  PDE -based  tractography  was  developed  and  presented  at  the  2012  IEEE  Workshop  on 
Mathematical  Methods  in  Biomedical  Image  Analysis  (1)  and  the  2012  Scientific  Meeting  of  the  International  Society  for  Magnetic 
Resonance  in  Medicine  (2).  Details  of  the  theory  and  fast  implementation  are  included  in  the  attached  publication  (1).  Highlights 
include: 

•  For  the  first  time,  encoding  of  termination  constraints  on  the  PDE  tractography  within  the  formalism  itself  as  opposed  to  ad- 
hoc  filtering. 

•  Acceleration  over  standard  probabilistic  tracking  by  several  orders  of  magnitude  from,  for  example  10  cpu-hours  for  a  single 
track  with  standard  approaches  to  2  cpu -seconds  with  the  new  approach. 

•  The  acceleration  of  the  implementation  over  standard  approaches  will  enable  whole-brain  tractography. 

An  important  feature  of  the  new  approach  is  that  it  successfully  delineates  tracks  when  standard  tractography  fails  (figure  1). 


Validate  noninvasive  measurements  of  connectivity  by  comparison  to  gold  standard,  invasive  electrophysiology  measurements. 
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A  total  of  5  subjects’s  imaging  and  electrophysiologic  data  were  examined.  Although  we  originally  proposed  examining  data  from  10 
subjects,  electrophysiologic  data  were  incomplete  in  all  but  5. 

The  overall  strategy  is  to  compare  measures  of  anatomical  connectivity  from  tractography  with  gold  standard  measures  of  connectivity 
from  electrophysiology.  Figure  2  indicates  the  location  of  subdural  electrodes  in  one  subject.  Connectivity  is  measured  by  first 
identifying  Broca’s  area  on  the  grid  with  a  speech  arrest  task  (electrodes  1 1  and  12  on  figure  2).  Cortico-cortico  evoked  potentials 
(CCEP)  are  measured  between  these  and  all  other  electrodes  to  identify  points  of  highest  connectivity. 

The  new  PDE-based  methodology  was  able  to  identify  connections  between  any  given  pair  of  electrodes  in  all  subjects  examined 
(figure  3).  A  sagittal  view  of  the  brain  with  overlaid  tracts  derived  from  the  fast  PDE  methodology  is  shown  in  figure  4.  The  seed  point 
for  all  tracts  is  in  the  Broca’s  region  of  the  left  inferior  frontal  gyrus.  The  target  points  are  electrode  contact  points  of  subdural  grids 
overlaid  on  the  left  parietal  and  temporal  lobes.  As  emphasized  by  the  color  coding,  the  tracts  can  be  grouped  into  three  bundles,  which 
correspond  to  known  physiology:  the  red  tracts  show  the  cluster  connecting  Broca’s  region  to  the  presumed  Wernike’s  region  of  the 
left  temporal-parietal  lobes  along  the  superior  longitudinal  fasciculus,  particularly  the  arcuate  fasciculus.  The  blue  tracts  fan  out  across 
the  external  capsule  to  connect  to  the  superior  and  middle  temporal  gyrus.  The  green  tracts  are  bundled  along  portions  of  the  left 
uncinate  fasciculus. 


Figure  2.  Photograph  of  subdural  electrodes  (left)  and  functional  maps  derived  from  electrical  stimulation  (right). 
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Figure  3.  Tracks  from  8  selected  electrodes  (red  circles)  to  all  other  electrodes  in  subdural  grid  (terminations  of  yellow 
lines). 


Figure  4.  Sagittal  view  of  the  brain  with  overlaid  tracts  derived  from  the  fast  PDE  methodology 


Correlation  with  electrode  measures  of  connectivity,  however,  was  modest.  Figure  5  shows  a  two  dimensional  scatter-plot  comparing 
electrophysiological  (EP)  connectivity  vs  DWI  connectivity,  in  an  epilepsy  patient  with  left-sided  subdural  grids.  Paired  electrode 
contacts  overlaying  Broca’s  region  in  the  left  inferior  frontal  gyrus  were  stimulated  using  a  CCEP  protocol  (8  mA,  1  Flz,  alternating 
unipolar  pulse  with  0.3  msec  duration).  Recordings  were  made  over  92  grid  contacts  placed  over  the  left  parietal  and  temporal  lobes. 
Tractography  was  performed  using  the  fast  PDE  approach  between  the  Broca’s  contact  and  all  remaing  electrodes,  for  a  total  of  92 
tracts.  The  connectivity  of  each  track  was  scored  in  a  measure  reflecting  tract-density,  and  these  values  are  scored  along  the  x-axis. 
The  y-axis  shows  the  magnitude  of  the  electrophysiological  response  between  10-20  msec  after  the  stimulus.  A  modest  correlation  is 
seen  (r=0.45),  many  contacts  show  either  a  strong  EP  connectivity  with  a  poor  DWI  connectivity,  or  a  poor  EP  connectivity  with  a 
strong  DWI  connectivity.  Furthermore,  the  correlation  does  not  achieve  statistical  significance  (p  >  0.2). 
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Figure  5.  Correlation  between  Tractography-based  and 
electrophysiologic  connectivity. 


This  result,  presented  at  the  2012  scientific  meeting  of  the  American  Society  of  Neuroradiology  (3),  indicates  that  the  link  between 
electrophysiologic  and  anatomical  connectivity  is  not  strong.  Competing  hypotheses  for  the  weak  correlation  include: 


1.  Assuming  that  tractography  identifies  true  axonal  connections,  electrophysiologic  connections  may  proceed  over  other,  less 
direct  routes. 

2.  The  weak  correlation  may  indicate  the  limits  to  reliability  of  tractography-based  delination  of  axonal  connections. 

3.  Coregistration  of  electrode  locations  to  image  positions  is  hampered  by  brain  shift  and  the  absence  of  post-implantation 
images. 


The  first  two  items  may  be  addressed  by  direct  examination  of  white  matter  pathways  on  pathology  from  surgically  resected  tissue.  In 
this  way,  it  may  be  possible  to  independently  test  if  tractography  identifies  white  matter  connections.  Also,  the  presence  of  an 
anatomical  connection  does  not  ensure  the  presence  of  electrical  activity  along  that  pathway.  Future  work  will  involve  the  use  of 
functional  connectivity  measurements  (4,  5)  to  incorporate  measurements  of  activity.  The  coregistration  issue  may  be  addressed  by 
taking  advantage  of  newly  developed  stereotactically  implanted  electrodes  (figure  6)  which  obviate  brain  shift  and  in  which  electrode 
location  can  be  known  to  the  limit  of  the  accuracy  of  the  stereotactic  frame. 
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Figure  6. 
Stereotactically 
implanted 
electrodes  (A)  do 
not  lead  to  brain 
shift  found  in 
subdural  electrodes 
(figure  6)  because 
the  skull  is  left 
intact.  Location  of 
electrodes  and  also 
be  determined  to 
high  precision  after 
implanation  (B,C). 


Summarize  and  publish  results. 

Work  described  here  has  been  presented  at  3  international  meetings  (1-3).  Manuscripts  are  currently  in  preparation  describing  the 
theory  and  implementation  of  the  tractography  and  of  correlation  with  electrode  recordings. 

KEY  RESEARCH  ACCOMPLISHMENTS 

•  Development  of  a  theoretical  basis  for  PDE-based  tractography. 

•  Implementation  of  the  tractography  that  is  fast  enough  to  enable  whole-brain  tractography. 

•  Finding  of  weak  correlation  between  gold-standard  electrode  recordings  and  anatomical  connectivity  based  on  tractography. 

REPORTABLE  OUTCOMES 

•  Presentation  at  2012  IEEE  Workshop  on  Mathematical  Methods  in  Biomedical  Image  Analysis  (1). 

•  Presentation  at  the  2012  Scientific  Meeting  of  the  International  Society  for  Magnetic  Resonance  in  Medicine  (2).. 

•  Presentation  at  2012  scientific  meeting  of  the  American  Society  of  Neuroradiology  (3) 

CONCLUSION 

We  have  implemented  a  methodology  for  tractography  that  is  fast  and  able  to  assess  connections  throughout  the  entire  brain.  This 
methodology  promises  to  enable  whole-brain  assessment  of  networks  of  brain  connections.  The  correspondence  between  these 
connections  and  electrophysiologic  activity  is  not  simply  one-to-one.  This  limitation  may  be  addressed  by  including  other  imaging 
measures,  such  as  functional  connectivity,  and  other  ground-truth  measurements  such  as  stereotactic  recordings  and  myelin  stains  of 
ex-vivo  tissue. 

As  a  scientific  or  medical  product,  the  work  accomplished  represents  a  step  towards  totally  non-invasive  evaluation  of  the  brain  at  risk 
for  epilepsy.  Such  an  evaluation  would  enable  rapid  evaluation  of  pharmacologic  interventions  and  development  of  new  therapies. 
Unfortunately,  although  victims  of  traumatic  brain  injury  are  at  high  risk  for  developing  epilepsy,  there  is  no  clear-cut  way  to  predict 
or  evaluate  strategies  for  treating  epileptic  seizures.  On  the  near  term,  surgical  intervention  for  phamacoresistant  epilepsy  often  relies 
on  highly  invasive  electrode  monitoring  that  is  an  option  for  only  the  most  highly  motivated  patients.  Progress  toward  noninvasive 
detection  of  targets  for  surgical  resection  would  relieve  the  burden  of  suffering  among  these  patients  while  opening  up  new  treatment 
options. 
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Abstract 

MRI  tractography  is  the  only  method  that  noninvasively 
maps  neural  connections  in  the  brain.  Interest  in  its  use 
for  diagnosis  and  treatment  of  neurological  disease  is 
growing  rapidly.  Probabilistic  tractography  provides 
quantitative  measures  that  can  be  interpreted  as  the 
strength  or  reliability  of  connections,  but  Monte  Carlo 
implementations  can  require  impractical  computation 
times  and  have  difficulty  identifying  connections  between 
distal  regions.  Here,  we  develop  a  generic  logical 
framework  for  probabilistic  tractography  with  minimal 
assumptions  that  lends  itself  to  solution  by  standard  finite- 
difference  methods.  We  demonstrate  an  implementation 
that  outperforms  Monte  Carlo  approaches  in  terms  of 
computation  time  and  identifying  distal  connections.  The 
generality  of  the  logic  and  the  speed  of  the  implementation 
indicate  the  potential  of  this  approach  for  real-time 
mapping  of  whole-brain  neural  connections. 

1.  Introduction 

Accurate,  noninvasive  maps  of  neural  connections  in  the 
brain  have  the  enormous  potential  to  improve  diagnosis 
and  treatment  of  neurological  disease.  For  example, 
neurosurgical  treatment  for  medically  refractory  epilepsy 
requires  identification  of  regions  of  abnormal  electrical 
activity  in  the  brain  for  resection.  Implanted  electrodes 
used  for  this  purpose  are  limited  in  size  and  can  therefore 
miss  the  target  while  imposing  risk  from  their  invasive 
nature.  A  compelling  rationale  for  tractography  is  to 
augment  or  even  replace  such  invasive  procedures.  By 
providing  noninvasive,  whole-brain  maps  of  neural 
connections,  tractography  may  be  able  to  identify  targets 
for  resection  that  the  electrodes  miss. 

Unfortunately,  tractography  in  its  current  state  cannot 
provide  reliable  whole-brain  maps  for  practical  use. 
Standard  deterministic  tractography  [1-3]  misses  many 
known  connections  [4].  Probabilistic  tractography 
identifies  additional  connections  [5]  but  can  require 
computational  times  that  are  too  long  to  impact  treatment. 


More  important,  tractography  can  generate  false 
connections.  Statistical  methods,  based  on  probabilistic 
tractography  for  filtering  these  false  connections,  are 
promising  [6,  7]  but  incur  additional  computational 
burden. 

In  this  paper,  we  develop  a  generic  logical  framework 
for  widespread  tractography  that  lends  itself  to  a  fast 
numerical  implementation.  By  taking  into  account  a 
minimal  set  of  logical  conditions,  it  is  possible  to  use  the 
rules  of  probability  theory  to  construct  a  solution  for  the 
density  of  tracks  in  every  voxel.  The  solution  lends  itself 
to  fast  numerical  solution  by  finite  difference  methods, 
such  as  simultaneous  over-relaxation  or  conjugate  gradient 
algorithms  [8].  We  show  such  implementations  offer 
orders-of-magnitude  speed  advantage  over  more  standard 
Monte  Carlo  implementations  of  probabilistic 
tractography,  particularly  in  geometries  challenging  to  the 
latter. 

2.  Theory 

Our  objective  is  to  count  the  number  of  fiber  tracks  in 
each  voxel  subject  to  the  conditions  that  all  tracks  start  in  a 
given  seed  voxel  and  all  terminate  in  a  given  target  voxel, 
without  crossing  given  boundaries.  For  example,  the 
superior  part  of  the  corticospinal  tract  can  be  identified  by 
generating  fiber  tracks  seeded  in  the  posterior  limb  of  the 
internal  capsule  and  terminating  in  the  motor  cortex.  Note 
that  “tracts”  refer  to  anatomical  nerve  pathways  while 
“trac/cs”  refer  to  the  mathematical  curves  generated  by 
tractography  algorithms. 

We  begin  by  defining  a  quantity  that  will  be  shown  to  be 
equivalent  to  the  desired  track  count: 

<p(j)  =  "the  number  of  tracks  originating  in  voxel  j 
subject  to  the  condition  that  the  tracks  reach  the  target 
voxel  before  hitting  either  the  seed  voxel  or  boundary 
region.”  Note  that  <p(j)  forms  a  subset  of  all  tracks 
leavingy;  that  is,  of  all  tracks  originating  from  j,  many 
hit  the  boundary  or  return  to  the  seed  before  hitting 
the  target. 
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The  originating  voxel  may  be  set  as  the  seed  voxel. 
Furthermore,  sets  of  seed  or  target  voxels  can  be 
considered  instead  of  individual  voxels.  For  the  sake  of 
brevity,  the  condition  will  be  referred  to  as  “reaching  the 
target.” 

The  assumption  that  tracks  are  conserved  leads  to  the 
solution  for  <p(j).  The  assumption  simply  states  that  tracks 
neither  appear  nor  disappear  at  any  points  but  the  seed, 
target  or  boundary.  A  consequence  of  this  assumption  is 
that  the  value  of  <p(j)  is  directly  related  to  the  value  of  <p(i) 
at  all  neighbor  voxels,  i: 

<P(j)  =  'Lp(i,j)p(i)  (1) 

i 


where  p(ij)  is  the  probability  that  a  track  moves  from  i  to 
j,  subject  to  the  condition  that  it  reaches  the  target.  If  we 
assume  that  the  tractography  is  a  Markov  process  (has  no 
memory),  <p(j)  is  not  only  the  number  of  tracks  that 
originate  in  j,  subject  to  the  conditions,  but  is  also  the 
number  of  tracks  in  j,  the  desired  quantity.  Solving  for  <p(j) 
will  therefore  achieve  our  objective.  Achieving  the 
solution  has  two  components:  the  logical  task  of  deriving 
the  probabilities,  p(ij),  and  the  numerical  task  of 
calculating  the  values  of  <p(j )  once  the  p( i,j)  have  been 
derived. 

The  logical  task  of  relating  p(ij)  hinges  on  its  definition 
as  a  conditional  probability.  The  condition  that  tracks 
reach  the  target  can  be  encoded  directly  into  the 
probabilities: 


P(i,j)  = 


p{Ar\B ) 
P(B) 


(2) 


where  B  is  the  proposition,  “track  reaches  the  target”  and  A 
is  the  proposition,  “track  moves  from  i  to  j .”  We  define 
k(ij)  as  the  probability  of  the  track  moving  from  i  to  j  and 
r(j)  as  the  probability  that  a  track  originating  in  j  reaches 
the  target.  Assuming  the  independence  of  X(i,j)  and  r(j), 
the  joint  probability  of  the  track  moving  from  i  to  j  and 
then  proceeding  to  the  target  is: 


p(Ar>B)  =  X(i,  j)r(j)  (3) 


while  the  overall  probability  of  reaching  the  target  from 
the  originating  voxel,  i,  is  r(i).  We  therefore  find: 

(4) 

r(i) 

Figure  1  diagrams  a  concrete  example  illustrating  the 
intuition  behind  equation  4. 

Although  the  probability,  r(j),  is  distinct  from  the 
previously  defined  quantity,  <p(J),  which  is  the  number  of 


Figure  1:  Cartoon  example  showing  application  of  conditional 
probabilities.  Consider  1000  tracts  originating  from  voxel  i.  Of 
these,  only  a  fraction  (40%  in  this  example,  or  r(i)  =  0.4)  reach 
the  target  voxel  (for  example,  track  “C”),  with  many  tracts 
hitting  either  the  boundary  (“B”)  or  returning  to  the  seed  (“A”). 
In  this  2-D  example,  there  is  a  25%  probability  for  each  tract  to 
step  into  one  of  the  four  orthogonal  directions  (/.(/,  j)  =  0.25  is 
constant).  Consider  those  tracts  that  step  into  the  adjacent  voxel 
j.  Of  those,  200  reach  the  target  (i.e.  r(J)  =  200/250  =  0.8).  We 
can  thus  see  that  the  probability  for  a  successful  track  to  reach 
the  target  from  i  via  j  is  p(i,  j )  =  200/400  =  (1000  x  0.25  x 
0.8)/(1000  x  0.4)  =  Mi,  J)  x  r(j)/r(i). 

tracks  in  j,  we  can  use  a  similar  conservation  argument  to 
construct  a  solution: 

r(0  =  '2_JMi,j)r(j)  (5) 

j 

The  probability  that  a  track  originating  in  i  reaches  the 
target  is  the  sum  of  the  probabilities  that  the  track  moves 
to  neighbors  j,  each  multiplied  by  the  probability  of  tracks 
originating  in  j  reaching  the  target.  This  simply  uses  the 
fact  that  if  i  is  not  the  target,  any  track  must  pass  through 
some  neighbor  before  reaching  the  target. 

As  k(ij)  indicates  the  likelihood  that  a  track  moves  from 
i  to  j  regardless  of  any  other  conditions,  it  is  an  entirely 
local  property.  This  quantity  can  be  estimated  from  the 
data.  For  example,  previous  work  has  derived  a  similar 
quantity  by  integrating  a  function  of  the  diffusion  tensor 
over  a  solid  angle  around  a  line  connecting  neighboring 
voxels  [9].  Alternatively,  integration  can  be  implemented 
over  a  function  that  accounts  for  complex  fiber  geometries, 
such  as  the  diffusion  or  fiber  orientation  distribution 
function  (dODF  or  fODF,  respectively)  [10,  11]. 

The  logical  derivation  is  now  complete.  From  the  data, 
we  derive  the  local  probability,  1,  which,  through  equation 
5,  allows  us  to  derive  the  probability  r.  The  probabilities  p 
follow  through  equation  4,  leading  to  the  desired  quantity, 
<P- 
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the  seed,  r  is  set  to  0  at  the  seed.  The  number  of  tracks,  <p, 
can  be  set  to  an  arbitrary  value  at  the  seed  and  is  set  to  0  at 
the  boundary.  In  order  to  avoid  the  target  acting  as  a 
source  of  tracks,  <p  must  be  set  to  0  at  the  target.  For  the 
numerical  approach,  we  applied  both  a  simultaneous  over¬ 
relaxation  method  (e.g.  Gauss-Seidel)  and  a  biconjugate 
gradient  stabilized  method.  The  latter  recasts  each  finite 
difference  equation  into  a  massive,  but  sparse,  matrix 
equation. 

3.  Methods 

Solutions  were  implemented  in  1DL  (ITT  Visual 
Information  Systems,  Boulder,  CO)  and  MATLAB  (The 
Mathworks,  Natick,  MA).  Test  data  sets  were  used  to 
compare  this  solution  with  a  Monte  Carlo  probabilistic 
tractography  routine  [12].  The  latter  is  essentially  a 
weighted  random  walk,  in  which  the  probability  of 
stepping  in  a  given  direction  is  determined  from 
underlying  diffusion  imaging  data.  The  test  data  sets 
included  numerically  synthesized  phantoms,  a  hardware 
phantom  from  the  fibre  cup  competition  [13]  and  in  vivo 
data.  The  in  vivo  data  was  acquired  under  a  protocol 
approved  by  the  Cleveland  Clinic  Institutional  Review 
Board.  Whole-brain  high  angular  resolution  diffusion 
imaging  (2.5  mm  isotropic  voxels,  61  non-collinear 
diffusion-weighting  gradients  with  b=1000sec/mm2  and  7 
b=0  acquisitions)  was  acquired  on  a  Siemens  TIM  Trio 
(Siemens  Medical  Solutions,  Erlangen)  followed  by 
iterative  motion  correction  [14].  For  the  hardware  phantom 


Figure  3:  Corrected  track  densities  in  the  hardware  phantom 
using  the  Monte  Carlo  (left)  and  the  finite  difference  approaches 
(center),  which  required  60  hours  and  less  than  1  second, 
respectively.  A  fractional  anisotropy  map  (right)  is  shown  for 
reference. 


Figure  4:  Corrected  track  densities  calculated  by  the  Monte 
Carlo  (left)  and  finite-difference  (right)  approaches,  overlaid  on 
fractional  anisotropy  maps  for  anatomical  reference. 

and  in  vivo  data,  fODFs  were  calculated  in  each  voxel  by 
spherical  deconvolution  [10]  with  user-independent 
optimized  regularization  [15].  Local  hopping  probabilities, 
1,  were  calculated  by  integrating  over  the  solid  angle  of  a 
vector  connecting  each  voxel  with  its  26  neighbors.  In  the 
numerical  phantoms,  the  probabilities  were  set  directly. 
For  each  test  data  set,  the  solution  was  also  found 
assuming  isotropic  probabilities. 
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4.  Results 

A  configuration  in  which  a  seed  and  target  are  near  the 
boundary,  such  as  those  found  in  proximal  cortico-cortico 
U-fiber  connections,  is  particularly  difficult  for  the  Monte 
Carlo  approach,  but  not  for  the  finite-difference 
implementation.  Figure  2  shows  an  illustration  via  a 
numerical  phantom,  a  simple  2 -dimensional  region  with 
isotropic  probabilities.  Seed  and  target  voxels  were  placed 
near  the  boundary.  Track  counts  were  normalized  to  the 
value  set  at  the  seed,  yielding  a  normalized  track  density. 
While  the  Monte  Carlo  and  finite  difference  approaches 
yielded  track  densities  that  differed  by  no  more  than  2%, 
the  finite  difference  was  faster  by  four  orders  of 
magnitude. 

Figure  3  shows  results  from  the  fibre  cup  hardware 
phantom.  A  simple  correction  was  implemented  in  order  to 
account  for  chance  connections  associated  with 
probabilistic  tractography.  Track  densities  calculated 
assuming  isotropic  hop  probabilities,  X,  were  subtracted 
from  those  calculated  using  probabilities  derived  from  the 
fODF.  This  correction  has  been  developed  more  fully  in 
the  Monte  Carlo  framework  [6],  Although  the  Monte  Carlo 
approach  required  over  60  hours  and  the  finite  difference 
approach  required  less  than  1  second,  the  latter  resulted  in 
a  smoother  map. 

In  vivo  data  is  shown  in  Figure  4.  While  the  Monte 
Carlo  approach  required  10  hours  of  computation,  the 
finite  difference  approach  required  only  2  seconds.  The 
huge  difference  in  computation  time  results  partly  from  the 
difficulty  with  which  the  Monte  Carlo  approach  is  able  to 
determine  track  densities  at  regions  distal  to  the  seed. 
Much  iteration  is  required  by  the  Monte  Carlo  method, 
while  the  finite  difference  approach  converges  rapidly  and 
at  the  same  rate  at  proximal  and  distal  sites. 

5.  Discussion 

We  demonstrate  that  a  minimal  set  of  assumptions  and 
logical  considerations  can  lead  to  a  fast  solution  for 
tractography.  Most  previous  probabilistic  approaches  [16, 
17]  considered  only  the  local  probability,  X,  of  moving 
between  neighboring  voxels.  A  Bayesian  approach  by 
Friman  et  al.  [18]  considered  conditional  probabilities,  but 
with  more  specific  conditions  associated  with  the 
trajectory  of  each  track,  particularly  as  influenced  by 
boundaries 

A  number  of  issues  remain  to  be  resolved.  Optimal 
calculation  of  the  local  probability  may  depend  on  the 
choice  of  dODF  or  fODF,  of  which  there  are  many  types 
[19].  Other  numerical  solutions  may  yield  even  greater 
performance  benefits.  Also,  the  correction  approach  using 
subtraction  of  isotropic  track  densities  is  ad-hoc  and  not  as 
well-developed  as  that  derived  in  the  context  of  Monte 
Carlo  approaches  [6], 


6.  Conclusion 

The  proposed  method  yields  results  comparable  to  more 
standard  Monte  Carlo  approaches  to  probabilistic 
tractography.  Flowever,  the  enormous  boost  in  speed 
suggests  the  possibility  of  determining  whole-brain 
anatomical  connectivity  in  a  clinically  relevant  amount  of 
time. 
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